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PQ . Abstract 

>— ' ■ Telomere are nucleo-proteins located at the end terminal of eu- 

Q I karyotic chromosomes which shorten at each cell division. Recent ev- 



X) 



idences points to telomere as a cell clock mechanism which limits the 
total number of cell divisions (Hayflick limit) and decelerates cellular 
^7' proliferation. In a recent work, simulations of a discrete time stochas- 

tic telomere regulated model for cell division produced growth curves 
similar to those observed in mesenchymal human stem cells and quan- 
^ I titatively very similar to Gompertzian Growth. 

In this work we considered the continuous time version of the 
stochastic telomere regulated model for cultured cell division. Instead 
of simulations, we obtained the analytic expression for the expected 
mean cell population. The growth curve obtained is still very close to 
^^ ■ a Gompertzian growth, within the limit of 2% of the final size. More- 

00 I over, Gompertzian growth is retrieved as a limit case. We also evaluate 

^D ■ the mean telomere length, which can be checked against cell culture 

data. 

Gompertzian growth is observed in tumors, animals and tissue re- 
. , generation but its biological foundations are still unclear. This work 

jrt ' adds a possible intrinsic control mechanism for this growth model. 
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1 Introduction 

Growth and replicative senescence are very often connected. In cultured 
cells growth is not observed indefinitely, the division rate slows down and 
ultimately ceases [26]. A cell clock mechanism based on the shortening of 
the replica with respect to the DNA template has been proposed in 1970s 
and was called marginotomy of DNA [39]. There are significant evidences 
that telomere act as a molecular counting device that regulates the number 
of cell divisions and limits further division after a critic length is achieved 
|25| . Replicative senescence of cultured cells is mainly imputed to telom- 
ere shortening, which is a mitosis-dependent process pj. Age dependent 
telomere shortening is also observed in adult somatic cells |46j 

Telomerase synthesizes human hexameric repeats (TTAGGG) onto telom- 
eric ends. This is the main mechanism to overcome telomeric losses of each 
cell division [35]. Telomere elongation independent of telomerase activity 
has been described in yeasts ^5\ but it is not yet documented in animals. 
Telomerase is present in human embryonic tissues, is not detected in most 
adult tissues, but is up regulated or reactivated in most of human cancers 
[45j . Telomerase activity correlated significantly with the Ki-67-expression 
(a marker of proliferative activity) both in neoplasia [19[ H8[ |36] and non- 
cancer cells [m [32]. It is interesting to note that expression of Ki-67 is 
routinely used in clinical oncology to access tumor aggressiveness. 

A transgenic mouse model (E/i-myc) of Burkitt's lymphoma demon- 
strates that short telomere combined with telomerase inactivity suppress 
tumorigenesis and this phenomenon is also observed when apoptosis is in- 
activated, indicating that short telomere negatively affect cellular prolif- 
eration [15) . So, telomere driven replicative senescence may be an anti- 
oncogenic mechanism, preventing replicative errors on the DNA to accumu- 
late indefinitely [47] . 

A very simple discrete-time stochastic model of telomere regulated growth 
has been proposed in [42]. Computer simulation of that model produced 
growth curves strikingly similar to Gompertzian growth. Gompertz model 
was originally proposed as an actuarial curve [22] to model mortality of an 
aging population. In fact, this approach still has applications in current sur- 
vival analysis [71ll2[[M]. One century after its creation, the model surpassed 
its original realm and was proposed as a biological growth curve [54[ 153] . 
Since then, this function has been successfully used for modelling animal 
growth [29l [30[ H] and regeneration [4l[52]. Moreover, since the pioneering 
work of Laird [31], Gompertz curve has been successfully used to model 
tumor growth [371 [Ml [5Ql [9l [JOl [ID [5] . 



The astonishing fact is not only that Gompertz model fits successfully 
in many cases of biological growth, but that it's biological foundation is still 
not fully understood. Several attempts were made in order to derive the 
biological explanation of Gompertzian growth model. Some efforts came 
from the use of allometric relationships and power-law equation [44], cell 
kinetics that employs quiescence of cells \23\ [27] , or maturity velocities over 
the cell population [2T] and entropy-based mechanism with thermodynamic 
characteristics [34J . 

Here, we take a step further on [42], by studying a continuous-time 
stochastic telomere-regulated growth model. Now, instead of simulations, 
we solved analytically the model and obtained the mean growth curve, which 
is strikingly similar to the Gompertzian growth. Moreover, Gompertzian 
growth is obtained as a limit of the stochastic model, when the mitosis 
number/Hayflick limit is high. We also estimated the mean telomere size of 
the cell population, which can be actually measured in cultured cells. 

Is worth of mention that the results obtained here do not depend on 
the particular division-counting mechanism. So, telomere or any other cell 
clock may be responsible for the replicative senescence/Gompertzian Growth 
modeled in this paper. 

2 The stochastic model 

We will assume that each cell in the initial cell population has its telomere 
with the same length Lq . For the sake of simplicity, we will also assume that 
a fixed amount of basis is lost by each telomere at each cell division, say 6. 
This is an approximation of the dynamics of telomere length [33]. Future 
research may consider more elaborated dynamics. So, after k divisions, the 
telomere length is 

L{k) = Lo- k6. (1) 

Our main biologic assumption is that, for each cell, division (mitosis) is an 
independent random process which probability depends only on the telomere 
length and duration of the time interval under consideration. Therefore, 
at any moment, time to division is a random variable with exponential 
distribution with parameter A which depends only on the telomere length 
L: 

A = A(L), P{t) = Ae^^* 

Probability of division shall decrease as the telomere shortens, and reaches 
for telomere length below a critical value L^i^. For the sake of simplicity. 



we will model this dependence as a linear one jl8]: 

\ _ J ( ~ -^minji L > Lmin 
[0, L < Lmin. 

Linear approximation is usual in natural sciences and indeed is a first order 
approximation, whenever we have differentiability of the modeled magni- 
tudes. 

To keep the model simple, we will assume that Lq — -Lmin is a multiple of 
5, which, in view of ([l]), means that when a cell reaches mitotic senescence, 
its telomere length is exactly Lmin- The maximum number of mitosis for 
the cell in the initial population is 

n = 



and 

Lq = Lmin + n5. 

So, using ([1]) and ([2]) we conclude that after k mitosis, the parameter A is 

Afc = 6(n — A;). 
Defining (3 as the parameter A for the initial telomere length Lq, we have 

/? = Ao = hn. 

Hence 

Afc = /3(n — k)/n, A; = 0, 1, • • • , n (3) 

Division is not an instantaneous process. Moreover, after a division each cell 
must synthesize a number of cellular components before it can divide again. 
Indeed, after a division, the cell mass is divided by 2. Even though, we will 
assume that these processes occur in a time scale much smaller that 1//3, 
which is the mean division time of the initial cell population. The model: 

• The initial state is a single cell with telomere length Lq. 

• This cell can undergo at most n mitosis. 

• Time to next mitosis for each cell is a random process with exponential 
distribution A^e" *=*, where k is the number of previous mitosis from 
the initial cell and 

Afc = /3(n — k)/n, A; = 0, 1, • ■ ■ , n — 1 



Let X{k,t) be the number of the cell population at time t which has 
undergone exactly k mitosis in his history. The expected value of X{k,t) is 



EiXik,t)) = Y,JPiX{k,t)=j). 

j=0 



Define 
Then 

and 



Xk{t) = E{X{k,t)), k = 0,---,n. 

Xfc(O) 



1, k = 0, 

0, k = 1,- ■ ■ ,n 



dxo{t)/dt = -XoXo{t), 

dxk{t)/dt = 2Xk^iXk-i{t) - \kXk{t), k = l,---,n. 

where A^ = 0. Note that 

Xk = (3{n — k)/n, /c = 0, 1, ■ ■ ■ , n. 

In matrix notation, omitting the zero entries, ((Tj) becomes 





Xq 






Xl 




d 






dt 


Xn-\ 
_Xn 





-Ao 

2Ao — Ai 

2Ai 



— A„-2 

2An-2 — A„_i 

2Ari„i — A„ 



Xo 
Xl 



^n— 1 



Hence, defining 



and M E i?("+i)x(«+i), 



x(t) = (xo(t),--- ,Xn{t)) 



M = {■mij}ij=i^...^n+i 



equations ©-([T]), becomes 



' '^'j 



-Aj_i, i = j, 
2Aj_i, i = j + l 
otherwise. 



(4) 

(5) 
(6) 

(7) 
(8) 



(9) 



x(0) = (l,0,--- ,0), ^x = Mx. 



(10) 



The solution of this hnear ODE is 



x(i) = exp(fM) 



(11) 



vahd for t > 0. The expected number of cells at time t is given by 

n 

J]x,(t) = [l,l,...,l]x(t). 



i=0 



3 The dynamics of the stochastic model 

In this section we present the analytical solution of the differential equations 
which governs the expected size of the cell population (jlOp . We compute in 
particular the matrix exp(tM). 

In order to solve (jlOp . first note that, since M is lower-triangular, the 
spectrum of M is — Aq, — Ai, • • • , —Xn- Define, for 7 G R and i,k = 1,2 .. . ,n+ 
1 

'0 i <k, 



a{i,k,-f) =a(7). 



i,k 



,1. l'n+1 - k^. . ^ , 



We claim that 



(12) 



(13) 



is a right-side eigenvector of M, corresponding to the eigenvalue — Afc_i. To 
check this fact we must evaluate 



(M + Afc_il)wfc 



where I is the {n + 1) x (n -|- 1) identity matrix. Using again the fact that 
M is lower-triangular and (J12p we have 



" bi,k ' 


b2,k 


pn+l,k_ 



bi,k = 0, i < k. 

For i = k, since mk,k = —Xk-i, 

bk,k = {'mk,k + Afe-i) Wk,k = 0. 



So, if k = n + 1 then 6 = and the claim holds. 
For the case k = 1,2,- ■ ■ ,n and i > k, 

bi,k = mi^i-i Wi^i^k + {mi^i + \k-i)wi^k 
Using also ([9]) and ([8]) we have 



\k 



2Ai_2 Wi^i^k + (-Aj-i + \k-i)wi^k 

(3 

-[2{n + 2 - i)wi^ik + {i- k)wik] 

n 



-2{n + 2-i) 
n 



i — k 

Wi-l,k + TTT -T, r^'Wi^k 

2[n + 2 — 1) 



Note that, according to (J12p for these indexes (i. e. A; < z < n + 1) 
Wi_i^k a{-2)i_i^k i - k 



Wi^k 



a(-2)i,fc -2(n + 2-i) 



The combination of the two above equations yields bi^k = 0, which concludes 
the proof of the claim on w^ being a eigenvector of M associated with the 
eigenvalue —Xk-i 

So, defining A(7),D G M{«+i)x{n+i)^ 



A(7) = {a(7)ij}, D = diag{-Ao,...,-An-i,-An}. (14) 



We have 
and 



M = A(-2)D [A(-2)] 



-1 



exp(tM) = A(-2) exp(iD) [A(-2)]-^ (15) 

To evaluate the above expression we will need two auxiliary results. 
Lemma 1 Let A(7) be as defined in (I12p . Then, for any r/ S M, 

[<,<~i,...,r/,l]A(7) = [(7? + 7)",(,? + 7)"-i,...,(r? + 7),l]. 

Proof. Let 

(m,--- ,Un+i) = [r/",<-\...,r/,l]A(7). 

Direct calculation, together with (fT2]l yields 

n+l 

Uk = ^V""^^'' ai^kh) 
1=1 

n+l 



E/^'-'^'-'i^'V-'k') 



To conclude the proof, define j = i — k and use the bmomial theorem to 
obtain 



n+l-fe 



Uk 



y^ iy"+l-*:-iy 



j=Q 



; ,Yn + 1 - /c 

J 



(r/ + 7) 



71+l-fc 



D 



Corollary 2 For any 7, /x G M, 

A(7)A(^)=A(7 + /i), 

and A(0) = I. In particular [A(7)]^ = A(— 7). 

Proof. The first equality follows trivially from Lemma[TJ The second equal- 
ity is obtained substituting 7 by in (J12p and using the definition of A (7) in 
p4p . The last part of the corollary follows directly from the previous ones. 

D 



Using Corollary [21 (fT5]) and ([TT]) we obtain 

x(t) = A(-2) exp(tD)A(2) 



(16) 



Note that e *''''= = (e *^/")" ^ . Therefore, using ([HD we have 

exp(tD) = diag{e"*^o, . . . , e~*^"} 

= diag{(e-*^/"r,---,(e-*'^/"),l}. 



(17) 



From this expressions and (jlip . we obtain a simple expression for the ex- 
pected number of cells at time t: 



Y, ^i{t) = [1,1,..., 1] A(-2) exp(tD) A(2) 



i=0 



(18) 



To evaluate (jlSp , use Lemma [J and (|17p to obtain 

[1,1,..., 1] A(-2) exp(tD) =[(-1)", {-ir-\ . . . , -1, , 1] exp(tD) 

= [(_g-/3t/")n^ (-e-^*/")"-\ . . . , -e"^*/", 1] 

Using again Lemma [J we obtain 

[(_e-/3t/n)"^ (_e-/3t/n)n-i, . . . , _e-/3t/n, i]A(2) = 

= [(2 - e-^*/")", (2 - e-/^*/")"-\ . . . , 2 - e'/^*/", 1] 
Hence, using the above equations and ([18]) we obtain 



4 = 

/ (,-l3t/n " 



2" 1 



Normalizing the final size to 1 we have 

5(t) = 2-"^x.(t)=(l-^^-j . (19) 

4 How close are the stochastic and Gompertzian 
models? 

First use (J19p to determine the time t* at which the (expected) population 
size is half of the final size: 

/ g-/3t*/n\ 

nln 1 =lnl/2 = -ln2 

Then 

/ e-/3t./n\ 
-Inh ^— j =(ln2)/n. 

As /i < — ln(l — /i) for < /i < 1, we have 

g-/3t./n 

< (ln2)/n. 



Hence, 

g-/3t*/n 



0„(ln2)/n, 0<en<l. 



2 
A Taylor expansion shows that 

n, f n \ log 2 , , 

where o„(l) vanishes as n | cx). 

Defining the normalized Growth curve 

^norm\T) ^ o (^t* -i- TIT) 

or, alternatively using the change of variable 

t = t* + nr, 
in ([19]) we get 

2-"f]x,(r) = (l - ^~^^^"^"^ y ^ exp (-e-^-e„ln2) . (20) 
i=o ^ n J 

Note that the last expression in the above equation is the Gompertzian 
growth curve. 

To verify how good is the approximation described in (j20p we analyzed 
an instance of the model with parameters in the range of those found in the 
somatic growth of animals and tumors. Human mesenchymal stem cells in 
culture can undergo around 30 mitosis before stopping division and usually 
do not express telomerase activity [55]. Here, we studied the model with 
n = 21. As the time scale is arbitrary, we set /? = 1. A Gompertzian Growth 
curve in this case, with final size 1 must have initial size 2~^^. So, consider 
the Gompertzian growth 

G{t) = exp(exp(-at)ln2"^^), 

where a is a parameter of the Gompertz curve, and is chosen to give a good 
fitness. So, we determined the best fitting parameter a = 1.116 and the 
maximum error for such q, (and 21 divisions) is 2%. Figure 1 compares 
both growth curves, using the best fitting a for the Gompertzian curve. 



10 




Figure 1: Stochastic and Gompertzian growth, 21 divisions 

5 Time evolution of the mean telomere length 

In cultured cells, telomere length is not evaluated individually. Instead, 
what is actually measured is the mean telomere length of a bunch of cells 
[31 [33]. Therefore, in order to verify the fitness of the stochastic model to 
real data, we must obtained the expected mean telomere length predicted 
by the model. 

In our model, the telomere length of a cell which has undergone k mitosis 
is 

Lmin + in- k)6. 

Therefore, the mean telomere length at time t is 

Ek=oiLmin + {n - k)6) X{k,t) _^ , , ELo(^-fc)^(fe>0 .^1^ 

Ek=0Xik,t) -^min + c) J^^^^^x{k,t) ■ ^''^ 

To evaluate this quotient, define for ^ : M x M — > IR 

n 

^{u,t) = ^u^-'X{k,t) (22) 

k=0 
11 



Let tp{u,t) be the expected value of ^(n, t) 

n 



(23) 



fc=0 



As E{X{k,t)) = Xk{t), using ([TH]) we obtain 



'ip{u,t) = [u^,u 



n „ n—1 



b",u"-i 



, c«. , , 



n,l]x(t) 

M,l]A(-2) exp(tD)A(2) 



Using Lemma [T] and ()17p we obtain 



r-u",?/"-! 



, lAJ , 



,u,l] A(-2) exp(tD) = 
[((n - 2)e-*^/")", ((u - 2)e-*'^/")"-\ • • • , (n - 2)6"*^/", 1] 

Combining the two above equations and using Lemma [1] again yields 



Note that 



V'('u,t) = (2 + (u-2)e-^*/"r. 

/ n \ n 

i? j;X(A;,t) =^Xi(t)=V(l,t). 



(24) 



^A;=0 



i=0 



and 



a 



^ X]^'^ " ^)^(^'*) = 5^(^ - ^)^iW = 7^V'(^,i)«=i 



a=o 



i=0 



Therefore, 



E Yl (^-- + (" - ^)'^) ^(^' *) =^mi„(2 - e-^*/")"+ 



^A:=0 



+ (5n(2-e~^*/")"~^e-^*/". 



The total telomere length of a cell population can be evaluated multiplying 
the mean telomere length of a sample by the population size. So, it it impor- 
tant to remark that this random variable can be measured in experiments 



12 



of cell culture tissue. Hence, we have an observable quantity which may 
agree or disagree with the model proposed here, confirming or invalidating 
respectively the model. 

Using the approximation 

ELo(" - k) X{k, t)\^E (ELo(^ - k) X{k, t)) 



^' n=oXik,t)) J' E{Y:U,X{Kt)) 
we obtain the estimation of the expected mean telomere size: 



E 



Y2=Q{L^n + {n-k)5)X{k,t) 



-t^min + 0' 



ne 



-I3t/n 



In Figure [2] we plotted the approximated mean telomere length as a func- 
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Figure 2: Mean Telomere Length as a function of mean population doubling, 
21 divisions, Lmin = 10.4kb, 5 = 86bp 

tion of the number of population doubles (log2 of the population) for some 
biologically feasible parameters [3]. This curve is similar to actual data of 
cultured mesenchymal stem cells presented in |3, Figure 3] 
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6 Discussion 

We studied an stochastic continuous-time growth model where the number 
of past divisions determine the probabihty of the next division. There is 
biological data supplying evidence that cell do count division, by means of 
telomere shortening. Moreover, mitosis-dependent telomere shortening is 
well documented, as well as the fact that cells with telomere length below a 
critical limit can not divide [16tl33j. Time lengths of S, G2 and M phases of 
the mitotic cycle are fairly constant. However, the duration of the cell cycle 
is variable from one cell to another in the same culture and most of this 
variability lies in Gl [2]. During Gl phase progression there is a checkpoint 
for cell senescence where telomere length is decisive [8]. Evidences that 
telomere shortening negatively affect the mitotic rate are still partial [15j . 
Here we supplied a theoretical evidence of this phenomenon, by proposing a 
model based on this assumption and which yielded the classical Gompertzian 
growth. Moreover, population doubling curves and mean telomere length on 
this model is consistent with actual data supplied by human mesenchymal 
cell culture [3] . Further validation of the model will be possible as soon as 
more accurate measures of mean telomere length of cultured cells become 
available. We regarded telomeric loss as a process with small variance and 
approximated it by a fixed shortening in each mitosis. 

In |43| telomere shortening was modeled as a random process with two 
different outcomes: a gradual telomere shortening or an abrupt shortening. 
The second outcome was related to sudden senescence. An aging mechanism 
based on the combination of telomere shortening, oxidative stress and DNA 
mutation was proposed in [l9]. Here, we modeled cell division as a stochastic 
process, which probability distribution depends on the telomere length. 

Modern cancer therapy is switching from the nonspecific toxic chemother- 
apy toward specific molecular targeted medicines. Malignant cell usually 
synthesizes telomerase, which prevents telomere shortening. If telomere 
shortening (down) regulates mitotic rate, this will have a significant impact 
in oncology, once specific telomerase inhibitors would act as very specific anti 
neoplastic agent [HI [T71 HQl [H]. For a general review on this perspective, 
see [28] 
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